Phase-field crystal modeling of equilibrium bcc-liquid interfaces 
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We investigate the equilibrium properties of bcc-liquid interfaces modeled with a continuum phase- 
field crystal (PFC) approach [K. R. Elder and M. Grant, Phys. Rev. E 70, 051605 (2004)]. A 
multiscale analysis of the PFC model is carried out which exploits the fact that the amplitudes of 
crystal density waves decay slowly into the liquid in the physically relevant limit where the freezing 
transition is weakly first order. This analysis yields a set of coupled equations for these amplitudes 
that is similar to the set of equations derived from Ginzburg-Landau (GL) theory [K.-A. Wu et ai, 
Phys. Rev. E 73, 094101 (2006)]. The two sets only difi'er in the details of higher order nonlinear 
couplings between diff'erent density waves, which is determined by the form of the nonlinearity 
assumed in the PFC model and by the ansatz that all polygons with the same number of sides have 
equal weight in GL theory. Despite these differences, for parameters (liquid structure factor and 
solid density wave amplitude) of Fe determined from molecular dynamic (MD) simulations, the PFC 
and GL amplitude equations yield very similar predictions for the overall magnitude and anisotropy 
of the interfacial free-energy and density wave profiles. These predictions are compared with MD 
simulations as well as numerical solutions of the PFC model. 

PACS numbers: 64.70.Dv, 68.08.-p, 81.16.Rf, 81.30.Fb 



I. INTRODUCTION 

The phase-field method is by now well-developed to 
simulate the continuum scale evolution of interfaces out- 
side of equilibrium with application to solidification 1] 
and other materials science problems [3, 3] . The method 
rests on a coarse-graining procedure that smears out the 
discrete atomic nature of the interface. Hence, the phe- 
nomenological form of the free-energy functional used 
to construct a conventional phase-field model generally 
needs to be tailored to reproduce quantitatively atom- 
istically determined interfacial properties. 

An important property in a crystal growth context is 
the anisotropy of the excess free-energy of the crystal- 
melt interface that is a key parameter controlling den- 
dritic evolution 0, H, [E S B B This anisotropy 
is traditionally incorporated phenomenologically into the 
phase-field model by letting the free-energy density de- 
pend on the direction normal to the interface, itself 
expressed in terms of the gradient of the phase-field 
[ill IT^ . fist . Computationally efficient implementations 
of these models have been successfully applied to simu- 
late dendritic evolution in materials with both atomically 
rough 0, H, [l3| and faceted [ij] interfaces. The conven- 
tional phase-field approach, however, falls short in prob- 
lems where crystalline defects have a profound influence 
on morphological evolution. For example, solidification 
twins can dramatically alter both eutectic [l^, [lB| and 
dendritic |l3 microstructures, and crystalline defects ul- 
timately control grain coalescence and microstructural 
evolution during and after the late stages of solidifica- 
tion 0. 



Over the last few years, the phase-field crystal (PFC) 
method has emerged as an attractive computational ap- 
proach to tackle this class of problems where atomic and 
continuum scales are tightly coupled [l^ [13, Hi], [12, HI] ■ 
This method is rooted in phenomenological continuum 
theories used to study equilibrium and noncquilibrium 
patterns with "crystal-like" ordering in diverse contexts. 
The models most closely related to the PFC model in 
their mathematical formulation have appeared in stud- 
ies of phase separation in block copolymers [IJj and 
Rayleigh-Benard convection [1^ [11] . 

Since the free-energy of the PFC model is a functional 
of the density of the material, the model can also be 
cast [m in the framework of classical density functional 
theory of freezing [27, 28, 29, 30, 31, 32, 33]. PFC simu- 
lations have the main advantage of resolving the atomic- 
scale density wave structure of a polycrystallinc material 
and of describing the defect-mediated evolution of this 
structure on time scales orders of magnitude longer than 
molecular dynamics (MD) simulations [l9l.l20Ll2lLl22l.l23l]. 

While the PFC method has been shown to describe 
qualitatively a wide range of phenomena [l^, HO, [ll|, 
m, HI], its predictive capability in a crystal growth 
context remains largely unexplored. We investigate in 
this paper to what degree the PFC model can repro- 
duce quantitatively some key equilibrium properties of 
the crystal-melt interface, in particular the magnitude 
and anisotropy of the interfacial free-energy 7. Well- 
developed atomistic methods to calculate these proper- 
ties [SJ, [H, m, [13, [H, [39I have been applied to both 
face-centered-cubic (fee) |38l . 
centered-cubic (bcc) systems 



and body- 
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Our study is based on the PFC model that is a re- 
formulation of the Swift-Hohenberg equation f2§\ with 
conserved dynamics introduced by Elder et al. (19, , ,20| . 
This model favors bcc crystal ordering in three dimen- 
sions. Our analysis of this model is closely related to pre- 
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vious studies of melting carried out in the framework of 
Ginzburg-Landau (GL) theory. The GL theory originally 
developed for bcc-liquid interfaces by Shih et al. 14p pre- 
dates the PFC model and was recently revisited j44"| in 
the light of recent results from MD simulations. This re- 
examination showed that GL theory yields predictions of 
7 and its anisotropy in reasonably good agreement with 
MD simulations for Fe, and the same MD results are used 
here to benchmark PFC model predictions. 

GL theory is derived from classical density functional 
theory that expresses the free-energy of the system as 
a functional of its density distribution n{r), as in the 
PFC model. Furthermore, it makes the strong assump- 
tion that n{f) can be expanded as a sum 

n(f) =no(^l + J2 "»(^^) ^'"^'"^ +•••)' (1) 

of density waves corresponding to the principal recipro- 
cal lattice vectors (where the index i spans the set of 12 
{-?^iio} vectors of the reciprocal fee lattice for bcc order- 
ing). The amplitudes Ui{r) of these density waves are the 
order parameters used to construct the GL free-energy. 
These amplitudes decay in the liquid at a rate that de- 
pends generally on the angle between Ki and the direc- 
tional normal to the solid-liquid interface, which makes 
7 anisotropic. The fact that the anisotropy predicted 
by GL theory is in reasonably good agreement with MD 
simulations suggests that this directional dependence is 
a main determinant of anisotropy [4^ . 

Since the crystal density field of the PFC model is 
also dominated by the principal reciprocal lattice vec- 
tors, we expect this model to yield similar predictions of 
bcc-liquid interfacial properties as GL theory. Of course, 
the two theories are not identical since the contribution of 
higher order reciprocal lattice vectors of magnitude larger 
than iXiiol, corresponding to " . . ." in Eq. ([1]) is small 
but non-vanishing in the PFC model. Furthermore, the 
strength of the nonlinear coupling between different den- 
sity waves is determined by the form of the free-energy 
functional in the PFC model, while it is determined in 
GL theory by using the simplifying assumption that all 
closed polygons composed of principal reciprocal lattice 
vectors with the same number of sides have equal weight 
[13,1411. Despite these differences, we find here that the 
PFC model and GL theory yield very similar predictions 
of bcc-liquid interfacial properties that are in reasonably 
good quantitative agreement with MD simulations. 

To relate formally the PFC model and GL theory, we 
carry out a weakly-nonlinear multiscale analysis of the 
PFC model. This type of analysis, pioneered in the con- 
text of Rayleigh-Benard convection [46] , has provided a 
fundamental understanding of the universal behavior of 
nonequilibrium patterns close to the onset of instability 
[2^ . It has also been revived recently in the framework of 
the renormalization group to derive computationally effi- 
cient implementations of the PFC model [13, HI] • In the 
pattern formation context where this analysis was first 



developed, the distance from the onset of instability can 
be characterized generally by a small parameter e, e.g. 
in Rayleigh-Benard convection e ^ {R — Rc) / Rc where R 
is the Rayleigh number and R^ is its critical value corre- 
sponding to the onset of instability. Furthermore, close 
to onset (e <C 1), spatially periodic patterns are generally 
slowly modulated in space. Considering the simplest case 
of a one-dimensional pattern for illustrative purposes, it 
is natural to write the field variables characterizing such a 
pattern in a form ~ A{Z)e^'^°'' + c.c.^ where Z ~ e^/'^z is a 
slow space variable, go is the wavenumber of the perfectly 
ordered pattern, and c.c. denotes the complex conjugate. 
The standard amplitude-equation approach consists of 
using a multiscale expansion to obtain an equation for 
the complex amplitude A[Z) starting from the under- 
lying equations governing the evolution of the pattern. 
The complex amplitude A{Z) = u{Z)e^'^'^^^ carries infor- 
mation about both the local real amplitude u[Z) of the 
pattern and its local spatial periodicity, or wavenumber 
q{Z) « go + e^^^^dz'^- Similarly, a dependence of the am- 
plitude on a slow time variable (omitted from the present 
discussion) can also be introduced to describe the slow 
temporal evolution of the pattern. 

For solid-liquid equilibrium, the pattern of interest is 
the three-dimensional crystal density field that is spa- 
tially modulated along the coordinate z normal to the 
solid-liquid interface. However, there is no direct ana- 
log of a small parameter e that can be made arbitrarily 
small by tuning some control parameter, such as the ex- 
ternally imposed temperature gradient in the example of 
Rayleigh-Benard convection. In contrast, e is uniquely 
determined by liquid structure factor properties when 
relating the PFC model to classical DFT. Thus e has 
a fixed value for a given material. For systems with low 
entropy of melting and atomically rough interfaces, how- 
ever, e turns out to be small enough (^ 0.1 for Fe) for 
a multiscale analysis to be just about justified quantita- 
tively. This smallness originates physically from the fact 
that density waves decay slowly in the liquid over several 
atomic layer spacings. This makes e, which is propor- 
tional to the square of the ratio of the layer spacing and 
the interface width, much smaller than unity. For faceted 
interfaces, however, density waves decay abruptly in the 
liquid and this expansion would break down. 

This paper is organized as follows. In section II, we 
briefly summarize the equations of the PFC model and 
construct the phase-diagram corresponding to bcc-liquid 
coexistence. In section III, we derive the amplitude equa- 
tions that describe the equilibrium profiles of density 
waves in the interface region from the aforementioned 
multiscale expansion. The phases $ of the complex am- 
plitudes turn out to be constant in the interface region at 
dominant order in this expansion, such that the density 
field can be described by Eq. (P) with real order param- 
eters that are the Ui{r)''s. This allows us to define the 
free-energy as a functional of these order parameters and 
to compare in section IV the PFC amplitude equations to 
GL theory [44| . This comparison is used to fix uniquely 
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FIG. 1: Phase diagram of the PFC model obtained under the 
approximation that the crystal density field is as a sum of 
density waves corresponding to the set of principal reciprocal 
lattice vectors for a given crystal structure. 



the parameters of the bare PFC model in terms of liq- 
uid structure factor properties and the solid density wave 
amplitude derived from MD simulations. Differences in 
the nonlinear coupling between density waves in the PFC 
amplitude equations and GL theory are also highlighted 
in this section. In section V, we compare quantitatively 
the predictions of 7 for different crystal faces obtained 
using (i) the direct numerical solution of the PFC model, 

(ii) the amplitude equations derived from the PFC model, 

(iii) GL theory [ij, and (iv) MD simulations. Finally, 
concluding remarks are given in section V. 



II. PHASE-FIELD CRYSTAL MODEL 



A. Basic equations and scaling 



uniquely fixed by the amplitude of density waves in the 
solid. 

To render the calculations less cumbersome, it is use- 
ful to rewrite the free-energy functional in dimcnsionless 
form by defining the parameter 



and making the substitutions, 



9_ 



F -^T, 



(3) 



(4) 



(5) 



(6) 



where all the transformed quantities to the right of the 
arrows are dimcnsionless and 



(7) 



In this study, we restrict our attention to equilibrium 
properties of the crystal melt interface. The condition 
that the chemical potential must be spatially uniform in 
equilibrium yields the equation 



5T_ 



(8) 



which is the starting point of the present study. Although 
the dimcnsionless formulation of the PFC model is more 
convenient to carry out calculations, we shall later trans- 
form the results back into dimensional form in order to 
make contact with GL theory and determine the phase- 
field parameters that appear in Eq. 



We consider the simplest PFC model defined by the 
free-energy functional [l^, [20| 



F 



(2) 



which is a transposition to crystalline solids of the Swift- 
Hohenberg model of pattern formation 25]. The con- 
served order parameter is a dimcnsionless measure of 
the crystal density field measured from some constant ref- 
erence value. The relationship of </> to the physical den- 
sity will be specified in the next section. The wavenum- 
ber qq sets the magnitude \Ki\ of the principal reciprocal 
lattice vectors that correspond to the first peak of the 
liquid structure factor S{K) at melting, and hence sets 
the scale of the ordered crystalline pattern ^ q^^. As 
shown in the next section, the parameters a and A can 
be related to properties of this peak and g, in turn, is 



B. Phase diagram 

To construct the phase diagram, we calculate sepa- 
rately the free-energy density (free-energy per unit vol- 
ume) as a function of the mean density in solid, de- 
noted by and hquid, //(^/i), using Eq. ([7]). We then 
use the standard common tangent construction, which is 
equivalent to equating the chemical potentials and grand 
potentials of the two phases, to obtain the equilibrium 
values of i/i in the solid (-0^) and liquid {ipi). 

Since the density is constant in the liquid, fi is ob- 
tained directly from Eq. ([7]) 



fi 



-(e-i)Tr + 



T 



(9) 



Furthermore, since e turns out to be a small parameter 
for spatially diffuse atomically rough interfaces, the solid 



4 



free-energy density can be well approximated by only 
considering the contribution of the principal reciprocal 
lattice vectors. Accordingly, the crystal density field can 
be written in the form analogous to Eq. ([T]) 



(10) 



4yls (cos qx cos qy + cos qx cos qz + cos qy cos qz) , 

where we have used the fact that all density waves have 
the same amplitude {\Ai\ = As) and all principal recip- 
rocal lattice vectors of the bcc structure have the same 
magnitude (|i?iio| = l-K^i-iol = •■• — V^q). The param- 
eters As and q are solved by substituting Eq. (fTU]) into 
Eq. ([7]) and minimizing the resulting free-energy with re- 
spect to As and q, which yields 



As 



2 7 1 , 



(11) 



and q — 1/ a/2, together with the expression for the solid 
free-energy density (for q = 1/ V2) 



fs 



-(6-1) 



IBip^Al - 



- 48V' + 135At (12) 



Applying the common tangent construction, which is de- 
tailed below in the small e limit, yields the bcc- liquid 
coexistence region in the phase diagram of Fig. [TJ Also 
shown are the other two-dimensional crystal structures 
(hexagonal and stripe phases) determined in previous 
studies using the same approximation where the crys- 
tal density field is a sum of density waves corresponding 
to the set of principal reciprocal lattice vectors [3, . 



III. DERIVATION OF THE AMPLITUDE 
EQUATIONS 

A. Small e analysis of the phase diagram 

For small e, we can seek a perturbative solution of the 
crystal density field of the form 



V'(r) = ^o(^) e'/' + ^i(^) e + e^^' + • 

and expand accordingly the average densities 

and 



(13) 
(14) 
(15) 



in the solid and liquid, respectively. Substituting these 
relations into the expressions for fs and fi, using the 
conditions of equality of the chemical potentials of the 



two phases, /^(V-s) =_f'i{i^i) =_IJ-e, 
grand potentials fs{ips) — t^'E'ips = 
collecting powers of e, we obtain 



^s 



and 



and equality of the 
fMi) - fJ'Etph and 



(16) 




(17) 



This shows that, in the small e limit, the PFC model 
exhibits a weak first-order freezing transition where the 
size of the solid-liquid coexistence region Atp = tps — ^'i ~ 
{ip^ — 'ipf)e^/'^ is much smaller than the mean value of 
the density ^ e^/^. These scalings imply that the mean 
density difference between the two phases only gives a 
small higher order correction to the density wave profiles 
through the interface and 7 in the small e limit. 



B. Multiscale expansion 



Using Eqs. p4)) and (|T5)) to evaluate the small e limit 
of the chemical potential ^lE — f'lii'i) = f'si'4's), the equi- 
librium equation of the density field ([S]) becomes 



(18) 



The derivation of the amplitude equation exploits the 
separation of scale between the width of the spatially 
diffuse interface and the interatomic layer spacing in the 
small e limit. This separation of scale allows us to assume 
that the envelope of density waves depends on a slow spa- 
tial variable Z = e^/^ z (i.e. Mr) = i^c + T, A°iZ)e'^'-^, 
and so on for higher order terms) where z denotes the 
coordinate along the direction normal to the solid-liquid 
interface. The multiscale expansion rests on treating the 
slow variable Z and the fast variable z as independent 
variables. Thus the spatial derivative along z transforms 
with the chain rule dz — > dz + e^^^dz, and the differential 
operator L"^ = (V^ -|- 1)^ in Eq. (fT5)) becomes 



^ + Ae^^^Ldz dz + 2e(L -I- 20^ 



(19) 



where the differential operator L on the right-hand-side 
only acts on the fast spatial variable z. 

Next, we substitute the small e expansion of the den- 
sity field into the equilibrium equation with the 
above transformation of the linear operator. Collecting 
terms with the same power e, we find at the order e^/^ 



which has the solution 

^^ = Y^A\[Z)t 



(20) 



(21) 
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where \Ki\ = 1 in our scaled units. At order e, we obtain 
lVi = 0, (22) 

which has the solution 

V^i=^Ai(Z)e^^^--, (23) 

i 

and collecting the terms at order e^/^ yields 

L^^2 + (4 dl - 1)^0 + V'o' - - V'c + V'c'. (24) 

The amplitude equations are obtained from the condi- 
tion for the existence of a solution of the above equation 
without needing to compute ■02 explicitly. Since L'^^p2 
gives a vanishing contribution for all density waves as- 
sociated with the set {A'^} of twelve principal recipro- 
cal lattice vectors of magnitude unity (i.e., L 



2AK-'. 



{-\K\' 

iiti 



I? 



iK-' 



— if \K\'^ = 1), all remaining terms 



e"'' ' must balance each other in order for a solution 
of Eq. p4)) to exist. For example, the condition that the 
coefficients of e'^on '- balance each other, yields 



(4(ifoii • nf dl + 3^,2 - l)Al,, + (3|AgnP + 6|A?io| 
+6|4lj2 + 6|^?oiP + 6|A?oil' + SKiH <i 

+"^011^101^101 + O^OlI^llO^lIO 



+6^ioi^iio'/'c 



= 0, 



(25) 



where everywhere in this paper z = fi corresponds to 
the direction normal of the interface that generally dif- 
fers from the crystal axes except for {100} crystal faces. 
This solvability condition must be satisfied independently 
for each Ki. This yields a set of twelve coupled ampli- 
tude equations (i.e., eleven additional equations to the 
one above) that are straightforward to obtain and we do 
not list them all here for brevity of presentation. These 
equations can also be obtained directly from the free- 
energy expressed as a functional of the amplitudes A^{Z) 
as described in the next subsection. 

The amplitude profiles are governed by these twelve 
coupled nonlinear amplitude equations. These equations 
can be reduced to a simple set of equations by consider- 
ing the symmetry of reciprocal lattice vectors. For {100} 
crystal faces, these twelve amplitudes can be separated 
into two subsets with the same value of {Ki-n)^ equals to 
1/2 and respectively. Therefore, the amplitude equa- 
tions are reduced to only two coupled equations and can 
be solved numerically. Similarly, we have two subsets of 
amplitudes for {111} crystal faces and three subsets of 
amplitudes for {110} crystal faces, which results in two 
and three coupled amplitude equations, respectively. 

As in GL theory 44] , the 7 anisotropy originates from 
the fact that the coefficients of the second derivative 
terms in the amplitude equations depend on (Ki ■ h) and 
hence on the orientation of the crystal face with respect 
to a fixed set of crystal axes. Furthermore, as mentioned 
in the introduction, since the amplitudes are complex. 



the spatial variation of the phase can cause the local wave 
vector to change through the solid-liquid interface by an 
amount proportional to the gradient of this phase. To 
determine this variation, we substitute 



^0(Z) = |A°(Z)|e^*-(^) 



(26) 



into the amplitude equations. We obtain that ^i(Z) = 
for the principal reciprocal lattice vectors that are or- 
thogonal to the interface normal, and 



1 



\Anz)\ 



dz{\AUZ)\H4Z))=0 



(27) 



for the other reciprocal lattice vectors. The above equa- 
tion implies that 



\A'dZ)\'dzMZ) - Co 



(28) 



where Co is a constant. Since the amplitudes must vanish 
in liquid, the divergence of d^i/dZ can be avoided only 
if Co ~ 0. Therefore, the wave vectors KiS are constant 
through the solid-liquid interface in the small e limit. 



C. Free-energy functional 

It is useful to express the free-energy of the solid-liquid 
system as a functional of the density wave amplitudes A^^. 
For this, we define AjF to be the free-energy measured 
from its constant value in the liquid. Since the ampli- 
tudes are non-conserved order parameters, the equilib- 
rium state simply corresponds to a minimum of this free- 
energy without extra constraint. This implies that AjF 
should be chosen such that the amplitude equations are 
recovered variationally from this free-energy. Namely, 
the equation for a given A^ derived in the last subsection 
should be equivalent to 



AT 



0, 



(29) 



up to a multiplicative constant. This constant can be 
determined by matching the limiting value of /V 
on the solid-side, where all the amplitudes are constant 
{A^ = e~^^'^As for all i) to the difference of free-energy 
densities between the two phases, fs — fi, where fs and 
// are given by Eqs. ^ and (fT^ and V is the volume. 
This yields the free-energy functional, 



AT = e^/^n / dZ 



dA° 



dZ 



+fiAr, 



(30) 



6 



where = J dxdy is the interface area and 



2^ 



EE 



'■\A^\ 



+0^110 ^110 ^101^101 

+ 0^110^011^011^110 

+0^011^101 ^101^011 
+6-0c^oii ^101^110 + 
+6-00^011 ^110^101 + 

i-DI/^c^Oll ^110^101 + 

+6V'c^oii ^?oi^iio + 



I fi^O AO AO * AO * 
+ "^110^110^101 ^101 

fi/lO * aO * aO * a" 
"^110 ^011 ^011 "^110 

I fi^o * AO AO * A^ 

6-00^011 ^?oi ^110 
fii/j A° A° 

OI/^c^Oll^llO ^101 
Oi/^c^Oll^llO ^101 

6i/'c^oii^?oi "^11 
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(31) 



It is simple to check that by applying Eq. ([29]) to the 
above functional for = ^oii we obtain the same am- 
plitude equation as Eq. and similarly for the other 
principal reciprocal lattice vectors. Finally, Eq. ^ im- 
plies that the dimensional free-energy functional derived 
from the amplitude equations (AE) is given by 



AFae 



AT. 



(32) 



IV. COMPARISON OF AMPLITUDE 
EQUATIONS AND GINZBURG-LANDAU 
THEORY 

In this section, we compare the free-energy functional 
derived from the PEC amplitude equations to GL theory 
p3 | . This comparison sheds light on the relation between 
this theory and the PFC model and uniquely fixes the 
parameters of the latter in terms of physical quantities 
that can be extracted from MD simulations. 



A. Ginzburg-Landau theory 

The free-energy functional of GL theory is expressed 
in terms of the amplitude ut of density waves defined by 
Eq. 11]). We write this functional here for convenience 
using the same notation as in Ref. [i^ l 



nokBTfl 



GL 



03 E/ '^y'' ^0,Ki+Kj+Kk 
-1-04 ^ CijklUiUjUkUi ^o^K.+K^+Kk+Ki 

J 2^ 

+b ^c. 



(33) 



du. 



dz 



where 5m,n is the Kronecker delta that equals or 1 for 
m n or m = n, respectively. The latter enforces that 
only combinations of principal reciprocal lattice vectors 



that form closed polygons Ki + Kj 



contribute 



to the free-energy functional. The multiplicative factors 
ai and b are introduced since it is convenient to nor- 
malize the sums of the c's to unity (i.e. '^^Ci — 1, 



'^ij^O,Ki+Kj ~ ^' etc). 

The coefficients of quadratic nonlinearities of the GL 
free-energy were determined in Ref. fi^l by relating 
AFqi^ to the free-energy functional that describes small 
density fluctuations of an inhomogeneous liquid in the 
simplest formulation of DFT (Eq. (3) in Ref. Ull) In 
particular, the latter can be reduced to the form [28[ 



AF, 



DFT 



dz 



E 



1 . 



dui 


2' 


dz 





(34) 



by assuming that the density wave amplitudes vary 
slowly through the interface region and are essentially 
constant on the scale of the interatomic layer spacing. 
The small e multiscale expansion of the last section is 
an alternative procedure to derive the form of Eq. ((34|) 
that formalizes this assumption. Here C (K) is the fourier 
transform of the direct correlation function C(|r|) 



C{K) = no / drC{\r\) 



(35) 



and S{K) = [1 — C(-ftr)]^^ is the liquid structure factor. 

Equating AFql and AFo ft at quadratic order in the 
nonlinearities and using the normalization that the sums 
of Ci 's and Cy 's equal unity, we obtain 



= 1/12, 



12 



02 



a 

~2C"{K„ 



(36) 
(37) 

(38) 
(39) 



where the magnitude \Ki\ — qq of the principal recipro- 
cal lattice vectors can be set equal to the K value corre- 
sponding to the first peak of the structure factor, K^ax, 
under the assumption that the wave vectors are constant 
in the interface region. This assumption was formally 
justified in the derivation of the amplitude equations in 
Sec. Ill B by showing that the phase $ of the complex 
amplitudes is constant in the interface region at leading 
order in the small e expansion. 

The reader is referred to Ref. [ij for the determination 
of the cubic and quartic nonlinearities in the GL theory, 
which shall be briefiy reviewed below. 
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B. Determination of phase-field crystal model 
parameters 

We are now in a position to compare the free-energy 
functionals derived from the ampUtude equations and GL 
theory and to relate the parameters of the PFC model to 
physical quantities. For this, we note that AFae has the 
same form as A.Fgl because the density wave amplitudes 
A'l's and Ui's proportionally related. The proportionality 
constant is readily obtained by combining Eq. ([5]) and 
Eq. ((131), which yields 



(40) 



Furthermore, Eq. ((29)) used to construct AFae is equiva- 
lent to the constraint that only combinations of principal 
reciprocal lattice vectors that form closed polygons con- 
tribute to the free-energy functional. 

Next, using Eq. ((40() and equating AFae, defined by 
Eqs. ((Sol) and ((Ml), and AFql defined by Eq. ((33]), we 
obtain the relations 




0-2 = 



12 no a(l 



12 



(41) 



C. Coefficients of quartic nonlinearities 



The free-energy functionals derived from the PFC am- 
plitude equations and GL theory only differ in the val- 
ues of the coefficients of higher order nonlinearities. As 
we shall see in the next section, these differences turn 
out to be unimportant because the amplitude equations 
and GL theory yield essentially identical predictions of 7 
and its anisotropy. However, they deserve brief mention. 
In GL theory, the coefficients 03 and a4 are determined 
from the two equilibrium conditions that (i) the solid and 
liquid phases must have equal free-energies at melting, 
AFGL\ui=u, = 0, and (ii) the equilibrium state of the 
solid is a minimum of free-energy, dAFg l I duj \u,=u.—^- 
These two conditions yields the relations [4J, Sal 



03 = 2a2lus 



(46) 



and 



a4 = a2 / u1 



(47) 



6=^^^ = -2C"(i^_.), (42) 

where we made use of Eqs. (1551) and ((M)) to write the 
second equalities and ■0c = — \/45/103 as shown earlier. 
Eqs. (dH) and ((i^ uniquely relate the parameters a and A 
of the PFC model to peak properties of the liquid struc- 
ture factor that can be computed from MD simulations 
or measured experimentally. They also fix the value of e 
related to a and A by Eq. ((3|) 



{l-i^c^)qlS{K^a.)C"{Kraax)' 

The only left unknown parameter g of the PFC model 
can be obtained by applying Eq. (|40p in the solid where 
all the density wave amplitudes have equal magnitude. 
Substituting into Eq. ((40)) Ui = Ug and the solid value of 
the amplitudes 

^0 ^ ,-1/2^^ ^ „^^^ + _Ly^rT^, (44) 

which follows from Eq. pip or Eq. ((^5)) . we obtain the 
relation 

.9 = A<z4(^-A^, + iy5^n^) /{nlu^J (45) 

This relation fixes g in terms of the other parameters 
and Us, which can be extracted directly from MD simu- 
lations [131 or related to the latent heat of melting and 
the temperature dependence of S{Kmax) 0,11^. 



The amplitude-equation free-energy functional AFae 
satisfies automatically the above two equilibrium condi- 
tions by construction. Thus, it only differs from AFql in 
the calculation of the other coefficients of the cubic and 
quartic terms, Cijk and Cijki ■ In GL theory, these coeffi- 
cients by the ansatz that all closed polygons of Ki's with 
the same number of sides have the same weight, which 
yields Cijk = 1/8 and Cijki = 1/27 

In contrast, in the PFC amplitude equations, these co- 
efficients are uniquely determined by the choice of the 
nonlinear terms in the original PFC free-energy func- 
tional. For the simplest choice of nonlinearity ~ cj)'^ con- 
sidered here, the amplitude equation derivation yielded 
the same coefficients of cubic terms as GL theory but dif- 
ferent coefficients of quartic terms. Comparing Eq. ()33p 
with Eqs. ()30p and ((5^ . we obtain that, in the expression 
for AF^^;, djki = 1/90 for two-sided polygons that con- 
tain only two wave vectors Ki and —Ki, and Cijki =4/90 
for the rest of the quartic terms. 

To make these differences explicit, we consider the 
{110} crystal faces. The set of 12 principal re- 
ciprocal lattice vectors corresponding to (110) direc- 
tion can be separated into three subsets with the 
same value of {Ki ■ fi)'^: subset I with 8 vectors 
([Oil], [Oil], [Oil], [101], [101], [101], [Oil], [101]) and {K, ■ 
hf = 1/4, subset II with 2 vectors ([110], [110]) and 
{ki ■ hf = 1, and subset III with 2 vectors ([110], [lIO]) 
and {Ki ■ nf — 0. Density waves in a given subset have 
the same amplitude denoted here by u, v, and w for sub- 
sets I, II and III, respectively. Then for {110} crystal 
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faces, Eqs. ((30|) and ([32)) reduce to 



AFae 



dz 



2 n 1 o In 

a,l-u +-V +-W 



1 2 1 2 ^ /^36 

-0-3 I 2" 2" V90"' 



16 



,2„,2 



16 



,2„,,2 



,,2„,2 



1 4,^ 1 . 

— w H It; 

90 90 

16 



H u f H u w H w V H u vw 

90 90 90 90 



-C"(|i?no|) 



du 
dz 



dv 


2' 


dz 





(48) 



and differ from the corresponding expression obtained 
from GL theory 



in section IIV C[ the density wave profiles were calculated 
by minimizing IS.Fae given by Eq. (|48p with respect to 
the order parameters u, v and w, and by solving numer- 
ically the resulting set of coupled ordinary differential 
equations with the boundary condition u = v ~ w — Ug 
in solid and u = v = w ^ Q in liquid. The value of 
7110 — I^Fae/^ was then computed by integration of 
Eq. ([^5)1 with these profiles. The same procedure was 
repeated for the {100} and {111} crystal faces, with dif- 
ferent set of order parameters for each crystal face. 

To compute 7 in PFC simulations, we first relax the 
density field -0 to a minimum of the free-energy functional 
T = J drf, where the free-energy density / is the inte- 
grand of Eq. ([7]), using a simple diffusive dynamics. We 
then compute 7 in dimensional units using the relation 



AFg, = ^^^^ / dz 



2 n 1 Q In 

a,\-u +-V +-W 



9 



dr 



27 



-w 



^11 ^22 "^22 ^2 

H u V H u w H w V H u vw 

27 27 27 27 



-C'di^nol) 



du 
dz 



C"{\Kno\) 



dv 


2' 


dz 





(49) 



V. COMPARISON OF CONTINUUM THEORIES 
AND MOLECULAR DYNAMICS SIMULATIONS 

In Ref. [44], the predictions of GL theory were com- 
pared to MD simulations of Fe with interatomic poten- 
tials developed by Mendelev, Han, Srolovitz, Ackland, 
Sun, and Asta (MH(SA)^) based on the embedded atom 
method [sot . In this section, we extend this comparison 
to include the predictions of both the PFC model, with 
the free-energy functional defined by Eq. ([2]) and the am- 
plitude equations derived from this model with the free- 
energy functional defined by Eqs. ([50)) and ([5^ . We use 
the same MD simulation results for the present compari- 
son. Details of the MD simulations and of the method to 
extract the density wave profiles from these simulations 
are given in Ref. [44j and need not be repeated here. 

The input parameters for the different continuum the- 
ories are computed from the MD simulations in order to 
make the comparison with these simulations as quanti- 
tative and precise as possible. These include the param- 
eters related to peak properties of the liquid structure 
factor K„,a^ = qo ^ 2.985 A~\ l/S{K^ax) = 0.332 , 

o 2 

C"{Kmax) = —10.40 A , and the amplitude of density 
waves corresponding to the principal reciprocal lattice 
vectors in the solid Us = 0.72. These input parameters 
fix the various coefficients of the continuum theories de- 
rived in the last section, which are listed in Table |T| 

The calculation of density wave profiles and 7 values 
for the PFC amplitude equations, Eqs. ([501) and ([5^ . 
proceeds in the same way as for GL theory [44| . For ex- 
ample, for the case of the {110} crystal faces elaborated 



(50) 

where i/is (V'/) and (/;) are the mean values of V' and 
the free-energy density in solid (liquid) , respectively, and 
fl = J dxdy is the interface area. Although e is small, 
these values need to be computed numerically (i.e. by 
calculating the solid free-energy density from the numer- 
ical solution of the PFC model rather than using the 
weakly nonlinear approximations derived in section II) 
in order to obtain an accurate computation of 7. 

The predictions of the different continuum theories are 
compared to MD simulations in Table II. Interestingly, 
despite the differences in quartic coefficients described 
in section IIV CI the PFC amplitude equations and GL 
theory give essentially identical predictions. The density 
wave profiles predicted by the two theories are almost 
indistinguishable on the scale of Fig. |21 Furthermore, 
the predicted 7 values by the different continuum theories 
for a given crystal face do not differ by more than a few 
tenth of a percent. Both theories predict a weak four- 
fold anisotropy £4 = (7100 - 7i 10)/ (7100 +7110) close to 
one percent consistent with the results of MD simulations 
with the MH(SA)2 EAM potential ^ for Fe. 

The PFC simulations predict essentially the same 
anisotropy value but about 10% larger 7 values that are 
in closer agreement with MD simulation results. The 
larger 7 values can be attributed to larger \K\ modes and 
to the variation of the mean density in the interface re- 
gion, both of which are neglected in the weakly-nonlinear 
amplitude equations and GL theory. 

The anisotropy parameter £4 defined in terms of 
7100 and 7110 has been traditionally used to quantify 
the magnitude of anisotropy in dendrite growth theory 
0, S H 0, 111- As seen in Table II, this parameter is 
reasonably well predicted by the PFC simulations and 
amplitude equations or GL theory. Over the past few 
years, however, numerous MD simulation studies have 
consistently found that at least two anisotropy parame- 
ters ei and 62 are necessary to represent the entire 7-plot 
of fcc-liquid and bcc-liquid interfaces in diverse systems 
Q. These parameters are defined by the expansion of 7 
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TABLE L Values of input parameters from MD simulations with interatomic EAM potential for Fe from MH(SA)'^ [H^l and 
resulting coefficients used in GL theory, the PFC model, and the amplitude equations derived from this model. 

no (A-^) aa b (A^) Us go (A'^) a (eV A"") X (eV A^) g (eV A^) e 
MD (MH(SA)^) (Ref. [41]) 0.0765 3.99 20.81 0.72 2.985 -2.136 0.291 9.705 0.0923 



TABLE II: Comparison of 7 values for different crystal faces (in erg/cm'^) and anisotropy parameters including £4 = (7100 — 
7iio)/(7ioo +7110) in percent and ei and €2 values (see text), predicted by MD simulations, and by various continuum theories 
(PFC simulations, PFC amplitude equations, and GL theory) with the input parameters of Table I from MD simulations. 





100 


110 


111 




ei 


£2 


MD (MH(SA)2) (Ref. [41]) 


177.0 (10.8) 


173.5 (10.6) 


173.4 (10.6) 


1.0(0.6) 


0.033 


0.0025 


PFC simulation 


160.47 


156.83 


152.00 


1.15 


0.075 


-0.0094 


Amplitude equations 


144.14 


140.67 


135.76 


1.22 


0.082 


-0.0110 


GL theory [44] 


144.26 


141.35 


137.57 


1.02 


0.066 


-0.0082 
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FIG. 2: Comparison of numerically calculated nonlinear order 
parameter profiles u and v for {110} crystal faces obtained 
from the PFC amplitude equations (dashed line) and the GL 
theory [3] (solid line) and computed form MD simulations 
with K\o\ and i?iio for u and v, respectively (solid circles). 



in terms of cubic harmonics (i.e., combination of spheri- 
cal harmonics with cubic symmetry) that has the form 



7(n) = 




^ ' (51) 



where the n^'s are the coordinates of the direction nor- 
mal to the interface (n) in a set of cartesian coordinates 
parallel to the crystal axes. Values of 7 for the the three 
independent crystal faces listed in Table II uniquely fix 
7o, ci, and £2. While a positive ei favors dendrite growth 
along the set of six (100) directions, a negative 62 favors 
growth along the set of twelve (110) directions. A recent 
phase-field simulation study has revealed the existence of 
hyper-branched dendrite morphologies with a basic set of 
twenty four growth directions between (100) and (110) 
over some region of the (£1,62) parameter space, where 
ei > and £2 < favor different growth directions (l0| . 



As seen from Table II, the agreement between the dif- 
ferent continuum theories and MD simulations is poorer 
for the ratio 7111/7100 than for 7110/7100- Consequently, 
the ei and £2 values, which depend on these two ratios, 
are not well predicted by these theories in comparison to 
£4, which depends only on 7100/7110- This discrepancy 
appears to be an intrinsic limitation of weakly nonlinear 
theories where anisotropy is computed using only one set 
of density waves associated with principal reciprocal lat- 
tice vectors of magnitude K^ax- While this one-set ap- 
proximation is reasonably good on the liquid side of the 
interface, where the density wave amplitude is small, it 
breaks down on the solid side where the highly nonlinear 
crystal density field is better approximated by sharply 
peaked Gaussians centered around atomic positions. Re- 
solving this field requires a very large number of sets of 
reciprocal lattice vectors [293. 

An interesting related issue is the sensitivity of crys- 
talline anisotropy to microscopic details of interatomic 
potentials. MD simulations to date indicate that the 
magnitude of this anisotropy tends to be larger for fee 
than bcc forming systems, suggesting that crystal struc- 
ture is a main determinant of anisotropy. Despite this 
trend, anisotropy values do depend on the choice of po- 
tentials for a given crystal structure. For example, two 
other interatomic potentials for Fe yield values of £4 twice 
smaller than for the MH(SA)2 potential [9]. 

In contrast, anisotropy values are independent of mate- 
rial parameters in both the PFC amplitude equations and 
GL theory. The reason is that all the material-dependent 
input parameters, which include the density wave ampli- 
tude in the solid and peak liquid structure factor prop- 
erties, can be scaled out of the free-energy functionals for 
these theories. This is readily seen in the dimensionless 
form of the free-energy functional for the PFC ampli- 
tude equations given by Eq. ([50)) . Consequently, the ra- 
tios of 7 values for different crystal faces that determine 
the anisotropy parameters £1 and £2 are universal for all 
bcc elements within the confines of each theory, and the 
value of anisotropy parameters depend on the nonlinear 
coupling between density waves. The results of Table 11 
show that differences in these couplings (i.e., coefficients 
of quartic terms in the free-energy functionals) lead to 
only small differences of anisotropy values. It is possible, 
however, that other choices could produce values of £1 
and £2 in closer agreement with MD simulations. 



VI. CONCLUSIONS 

We have studied equilibrium properties of bcc-liquid 
interfaces in a physically motivated small e limit of the 
PFC model [3, US] where the freezing transition is 
weakly first-order. This limit lends itself naturally to 
a multiscale analysis that was used to derive a set of 
equations for the leading order amplitudes of density 
waves corresponding to the set {Ki} of principal recip- 
rocal lattice vectors, and to express the free-energy of 
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the solid-liquid system as a functional of these ampli- 
tudes. Furthermore, by exploiting the close analogy be- 
tween this functional and GL theory derived from classi- 
cal DFT [S, '45*1 , we have determined all the parameters 
of the PFC model in terms of peak properties of the liquid 
structure factor and the solid density wave amplitude. 

In both the PFC amplitude equations and GL theory, 
the anisotropy of 7 originates from the directional de- 
pendence (i.e., the dependence on Ki ■ h where fi in the 
interface normal) of the coefficients of gradient-square 
terms (|VA°p) in the free-energy functional, which gov- 
ern the spatial decay rate of density waves in the liquid. 
In the isotropic limit where this directional dependence 
is neglected, and hence all amplitudes are equal, = (f) 
for all i, both theories reduce to the conventional phase- 
field model of solidification formulated in terms of the 
non-conserved order parameter (f>. From this standpoint, 
the present analysis relates formally the crystal and con- 
ventional phase-field models. 

Numerical results show that the PFC model, the am- 
plitude equations derived from this model, and GL the- 
ory [44] all give very similar predictions of 7 and its 
anisotropy for parameters of Fe where e is small enough 



(e « 0.1) for the amplitude equations to be quantita- 
tively valid. The magnitude of 7, the shape of density 
wave profiles in the spatially diffuse solid-liquid interface 
region, and the standard crystalline anisotropy parame- 
ter 64 defined in terms of the ratio 7110/7100 1 a-re in good 
overall agreement with the results of MD simulations. 
The various continuum theories, however, do not predict 
accurately higher order anisotropics that also depend on 
the ratio 7110/7100- These anisotropics probably depend 
generally on the contributions of higher sets of recipro- 
cal lattice vectors, which are neglected in the simplest 
formulation of the PFC model considered here. 
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